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We study structure formation in the presence of primordial non-Gaussianity of the local type with 
parameters /nl and (Jnl- We show that the distribution of dark- matter halos is naturally described 
by a multivariate bias scheme where the halo overdensity depends not only on the underlying 
matter density fluctuation 5 but also on the Gaussian part of the primordial gravitational potential 
ip. This corresponds to a non-local bias scheme in terms of 5 only. We derive the coefficients of the 
bias expansion as a function of the halo mass by applying the peak-background split to common 
parameterizations for the halo mass function in the non-Gaussian scenario. We then compute the 
' halo power spectrum and halo-matter cross spectrum in the framework of Eulerian perturbation 

^S) , theory up to third order. Comparing our results against N-body simulations, we find that our 

model accurately describes the numerical data for wavenumbers k < 0.1 — 0.3 h Mpc~^ depending 
' on redshift and halo mass. In our multivariate approach, perturbations in the halo counts trace 

on large scales and this explains why the halo and matter power spectra show different asymptotic 
trends for fc — >■ 0. This strongly scale-dependent bias originates from terms at leading order in our 
expansion. This is different from what happens using the standard univariate local bias where the 
scale-dependent terms come from badly behaved higher-order corrections. On the other hand, our 
biasing scheme reduces to the usual local bias on smaller scales where \'p\is typically much smaller 

■ than the density perturbations. We finally discuss the halo bispectrum in the context of multivariate 
' biasing and show that, due to its strong scale and shape dependence, it is a powerful tool for the 

• , detection of primordial non-Gaussianity from future galaxy surveys. 

Qh' PACS numbers: 98.65.Dx, 98.80.Cq 

o 

I. INTRODUCTION 

, Measurements of the cosmic microwave background (CMB) anisotropies have confirmed the hypothesis that the 
J — ' present inhomogeneities in the matter density were seeded by small fluctuations at primordial times Such per- 
turbations, which are expected to be created as quantum vacuum fluctuations, have been generally modeled with the 

' simple statistical assumption of being a Gaussian random fleld with nearly scale invariant power spectrum ^] . 

\ The inflationary mechanism is often invoked to describe the early universe, but the details remain debated 0]. 
^-H ■ In the simplest single-fleld, slow-roll model, small curvature (adiabatic) perturbations are generated with a nearly 
T-H ' Gaussian distribution 4, 5] . However in other models such as the curvaton scenario J6h8| some additional fields would 
On . decay at later times, producing larger non-Gaussianity [ol-fTTj: cyclic or ekpyrotic universes without infiation could 

■ also produce large non-Gaussianities during their contract ing phase flp. fisj . Furthermore, rnulti-field models can in 
^ , general produce isocurvature modes of the perturbations jl4| . See [15| for a review and [l6j for recent updates and 

future prospects. 

. ^ The first observable predictions from infiation — fiatness and the near scale invariance of the power spectrum of 

■ the perturbations — have been under scrutiny for some time from observations of both the large scale structure (LSS) 
" " ' [13] and the CMB f^. Fairly strict constraints on the adiabaticity of the primordial perturbations also exist jlSj, 

while their Gaussian distribution has onl y re cently become testable. 

Although other possibilities exist (see [l9| for a review), many models produce primordial non-Gaussianity of the 
local type where the Bardeen's potential <& can be expressed in terms of an auxiliary Gaussian potential Lp as 

oo 

<i>(x) = v'(x) + ^NL, [v" (x) - (X))] , (1) 

where the series will be in practice truncated at some finite order N , and the odd momenta of the Gaussian potential 
Lp vanish by definition. The first parameter Qnl2, usually dubbed /nl, quantifies the leading-order departure from 



X 



'Electronic address: IgiannantATastroDOTuni-boniiDOTde, porcianiATastroDOTum-bonnDOTde| 



2 



purely Gaussian initial conditions through the irreducible three-point function or the bispectrum of the potential. 
While standard inflation forecasts a slow-roll suppressed, primordial IfmJ ^ 1, subsequent evolutionary processes are 
expected to increase the amount of non-Gaussianity up to |/nl| ~ 1 [20|- On the other hand, more complex models 
can produce |/nl| 3> 1, although the actual predicted values vary. The second parameter Qnls, generally called ^nl, 
quantifies the next higher order contribution and is related to the irreducible four-point function or the trispectrum 
of the potential. Since (p ^ 10~^, this contribution can be important only if onl is big, ^nl ^ /nl- This is plausible 
in the interactive curvaton model [2T'-'23l| and in other multi-field scenarios 2J, 2^. 



The traditional method to constrain primordial non-Gaussianity has been the three-point statistics of CMB 
anisotropics. The current limits from WMAP are —9 < /nl < HI at the 95% c.l. !!!]. Different analyses of the 
same data found —178 < /nl < 64 using Minkowski functionals [ij, —4 < /nl < 80 using an optimized estimator 
[26j , and — 18 < /nl < 80 from wavelet decomposition [27|, while a detection (27 < /nl < 247) was claimed by [2^. 
Constraints on ^nl are —5.6 • 10^ < f/NL < 6.4 • 10^ [2^. The Planck satellite is expected to reduce the uncertainty to 
"■(/nl) ^ 5 [30| . This result will be nearly cosmic- variance limited, and further significant improvements from CMB 
studies will be difficult to achieve. This reason, together with the desire of having independent results, affected by 
different systematics, provided the motivation to study the detectability of primordial non-Gaussianity from the LSS. 
In this case, however, the non-linear growth of density perturbations can superimpose a new non-Gaussian signal onto 
the primordial one [311] , which may be difficult to retrieve. Determi ning the mass distribution of galaxy clusters at 
low and high redshift provides a way to circumvent this problem 's^, 'ssf! However, due to the low-number statistics, 
these methods have been so far less successful than the CMB. Upcoming surveys such as PanSTARRS, DES, LSST, 
ADEPT, EUCLID, JDEM or eROSITA, WFXT and SPT are expected to substantially improve the situation. 

A new technique, based on linear perturbation theory, has been recently introduced by [3J| (Dal07 in the follow- 
ing). These authors showed that local non-Gaussianity breaks the independence of small and large scales density 
fluctuations. As a consequence, the clustering of dark matter halos is altered, becoming enhanced on large scales 
for a positive /nl- An analytical derivation of the corresponding scale-dependent bias has been also presented by 
|35l - l38| . together with some observational constraints on /nl from existing data of the clustering of galaxies and 
their correlation with the CMB anisotropics. Using luminous red galaxies and quasars from the SDSS, Slosar et al. 
[36| obtained —29 < /nl < 69, competitive with the CMB results. The first constraints on ^nl from the LSS give 
-3.5 • 10^ < ffNL < 8.2 • 10^ ;39], assuming /nl = 0. 

N-body simulations show only approximate agreement with the model by Dal07 (40l - |4^ . In particular, the power 
spectrum of dark-matter halos seems to scale with the wavenumber and the /nl parameter in a different way than 
predicted. This discrepancy provides the main motivation for this paper where we study the effect of non-Gaussian 
initial conditions on the clustering of halos in the weakly non-linear regime of perturbation growth. Applying the peak- 
background split technique, we show that the halo distribution on large scales is naturally described by a bivariate 
local biasing scheme, where the halo overdensity is expanded in a Taylor series of both the matter perturbations S and 
the primordial Gaussian potential (p. Since ip and S are related by the Poisson equation, this can be equivalently seen 
as a non-local description in terms of S only. This reduces to the usual univariate local bias (where halo overdensities 
are expanded in series of S only) for Gaussian initial conditions and, in general, on small scales. Using standard 
(Eulerian) perturbation theory (SPT) up to third order to account for the non-linear growth of density fluctuations, 
we show that our new biasing scheme leads to the presence of several new terms in the halo power spectrum and 
bispectrum. We show that our results reduce to the usual Gaussian solution in the limit /nl — > 0, and to the results 
by Dal07 (revised as in 36, 40]) if we only consider the leading-order terms. Finally, we compare our theory with the 
N-body simulations by j4l| (hereafter PPH08) and find that our model can explain the numerical results to a much 
greater accuracy than both linear and univariate local theories. 

Our paper differs from the recent work by (43l - |46j which is based on the assumption that fluctuations in the 
halo number counts only depend on the local mass density. By ignoring the expansion in the potential ip and only 
considering the dependence on the matter density perturbations (5, one would obtain different results which do not 
reduce to the model by Dal07 to leading order and do not match the simulations as well; in this case higher-order 
terms in the halo power spectrum grow bigger than the first-order contribution on large scales, thus casting doubts on 
the validity of the perturbative expansion. Our approach is also different from the work by [47[ because we use SPT 
without applying any renormalization technique. Renormalizing the bias removes any undesired dependence on the 
cutoff scale introduced to regularize loop corrections in SPT; however, in such a model the bias coefficients cannot be 
predicted and should be used as fitting parameters to match observations or simulations. 

The plan of this paper is as follows. In Section|lT]we summarize the main biasing schemes which have been proposed 
to describe the distribution of different tracers of the LSS. In Section Hill we introduce some models for the halo mass 
function arising from non-Gaussian initial conditions and compare them against the N-body results by PPH08. We 
then describe in Section HVl how a multivariate bias scheme naturally emerges by applying the peak-background split 
technique to compute halo overdensities in the non-Gaussian case. In Section |V] we give a short summary of the 
statistical properties of non-Gaussian density fields. After computing the halo power spectrum and the halo-matter 
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cross spectrum in our multivariate biasing scheme in Section IVIl we test our theoretical models against the N-body 
simulations by PPH08. We derive the halo bispectrum in Section IVlIl and conclude in Section rVIIII 

II. TRACERS OF THE LARGE-SCALE STRUCTURE AND BIASING 

The large-scale structure of the Universe can be described in terms of different tracers: mass, luminosity, galaxy 
counts. In this paper we will consider dark-matter halos and mass but our formalism can be straightforwardly 
extended to any other tracer. Let us consider the mass overdensity field i5(x) and the corresponding density contrast 
of dark-matter halos in a given mass range (5;i(x). After smoothing both fields on a relatively large scale _R, it is 
reasonable to expect that Sh is a local function of S that can be expanded in a Taylor series as follows: 

Sh{^) = 60 + 6i(5(x) + 625'(x)/2! + 63<5^(x)/3! + ... , (2) 

where the bias coefficients bi are in principle scale and mass dependent [48l| . This approximation neglects stochasticity 
in the 6h vs. 6 relation and is thus dubbed local deterministic biasing. Numerical simulations from Gaussian initial 
conditions show that it is accurate on scales of the order of 10 Mpc and larger [i^. In the following sections we will 
show that Eq. ^ does not hold in the presence of primordial non-Gaussianity of the local type and we will explain 
how it should be modified. 

Note that, in general, the bias coefficients in Eq. ([2]) are not independent as the mean halo overdensity must vanish 
and Sh must assume the value —1 when 6 = —1. In order to build a predictive theory, the values for the bias coefficients 
should be derived from a model. A common approach is to use the peak-background split technique [1, [50l - l52j where 
the mass perturbations are divided into fine-grained (peak) and coarse-grained (background) components. The key 
idea is to ascribe halo formation to the collapse of the high-frequency modes, while the large-scale distribution and 
motion of these condensations is determined by the low-frequency modes. Starting from a model for the conditional 
halo mass function (i.e. the mass function in regions where the background density assumes a specific value), the 
peak-background split gives an expression for the halo distribution in Lagrangian space (i.e. in the linear density 
field, 61): 

5^(q) ^b^ + 6[<5i(q) + b^Sf{ci)/2l + b^6l{q)/3l + ... , (3) 

where the bias coefficients 6f are obtained from the i-th order derivatives of the conditional mass function with respect 
to the background density contrast (see Section HVl for further details). When the background scale is much larger 
than the Lagrangian size of the halos, the Lagrangian bias parameters show very little dependence on the background 
scale and the unconditional mass function can be safely used to derive them [53] . We will follow this approach in this 
paper. 

In the absence of large-scale velocity bias, the halo density in the evolved Eulerian space is given by 

l + <5,(x) = [l + <5^(q)][l + <5(x)] (4) 

where q is the Lagrangian position of the fluid elements that moved to the Eulerian location x [s^] ■ Note that the 
conversion between Lagrangian and Eulerian quantities is non-local, non-linear and stochastic as it depends on the 
displacement x — q and on both the initial and the evolved fields Si and S. Therefore, the local Lagrangian biasing 
scheme given in Eq. ([3]) will not generally be compatible with Eq. ([2]). Catelan et al. [H^ showed that these two 
bias models give rise to different shapes of the halo bispectrum that could then be used to distinguish between them 
using data from observations or simulations. A simplified approach is obtained by assuming that the long- wavelength 
modes of the density field evolve locally according to the spherical collapse model (sil. [ssj. In this case, Eqs. ^ 
and ^ are fully compatible and the Eulerian bias parameters can be written in terms of the Lagrangian ones (see 
Eq. in Section ITV)). In particular, bi — 1 + fef . This equation is completely general as it derives from mass 
conservation fsii, i53|. However, the relation between higher-order Eulerian and Lagrangian bias parameters depends 
on the adopted dynamics for the background density field. A perturbative calculation of the power spectrum for local 
Lagrangian biasing in the Gaussian scenario has been presented by j56l |. The equivalent result for the local Eulerian 
bias scheme has been derived by [H^l- In this paper we generalize this latter result to non-Gaussian initial conditions 
of the local type and also present a model for the bias coefficients as a function of the halo mass. The derivation of 
a multivariate bias scheme and the corresponding calculations of the halo power spectrum and bispectrum constitute 
our main results. 
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III. HALO MASS FUNCTION AND PRIMORDIAL NON-GAUSSIANITY 

The number density Af of halos of mass AI at a redshift z is described by the mass function 

d hi cr-i 



dln(A/) 



(5) 



where 6c — 1.686 is the threshold for the linear density contrast which corresponds to the collapse of spherical 
perturbations. In Eq. ([5]), denotes the variance of the linear density field, calculated as 



(M,z) 



27r2 



Po{k)Wf{k,M)dk, 



(6) 



with Po{k) the linear matter power spectrum at z = 0, D{z) the linear growth factor of density fluctuations normalized 
to unity today, and Wf{k,M) a filter function with mass resolution M. We use a top-hat filter in real space with 
radius Rf = [3M/(47rp)]^/'^, where p is the average density of the Universe. The analytical form of the distribution 
f{5c/<j) can be derived from a theoretical model or by fitting numerical data. Wc list below some possible choices for 
this function. 



Gaussian mass functions 



Press-Schechter (PS) For reference we first consider the Press-Schechter theory (58|, in which the mass function 
deriving from Gaussian initial conditions is given by 



/PS 

a 



2 (5c -A 
e 2ct- 

TT a 



(7) 



It is well known that this model gives only a rough approximation to numerical data (see Fig.[T]). The Press-Schechter 
theory can be improved by introducing extra parameters in the mass function and fitting them against numerical 
simulations. This approach has been followed e.g. by Jenkins et al. [so']. Warren et al. (W) [691, Tinker et al. i61i |. 

Sheth-Tormen (ST) The Press-Schechter theory is based on the spherical collapse model. This can be improved 
upon by considering the collapse of ellipsoidal perturbations and fitting some new parameters against numerical 
simulations, jssl. 1621. The final result, known as the ST mass function, is 




1+ 



5c - 
— e 

(7 



(8) 



where the extra parameters are a 
halos, which gives A = 0.322. 



0.707, p — 0.3 and A is obtained by requiring that all the mass is collapsed into 



B. Non- Gaussian mass functions 

In the simplest case, local non-Gaussianity is described by truncating Eq. ([T]) after the second order term, which 
corresponds to: 



i'(q) = ^(q) + /nl [^^(q) - (^2(q)>] 



(9) 



It can be shown that the halo mass function is very sensitive to the value of /nl- For positive (negative) values of 
/nl, its high- mass tail becomes more (less) prominent than in the Gaussian case. To first-order in the non-linearity 
parameter /nl, it is possible to account for this effect by considering the skewness of the density perturbations, defined 
as 53(a) - {5^)/(T^. 

Matarrese- Verde-Jimenez (MVJ) The Press-Schechter theory can be generalized to non-Gaussian initial con- 
ditions by using the saddle point approximation to calculate the probability for the linear density field to be above 5c 
1321. In this case, one obtains: 



/mvj — 

(T 



51 dSsja) ^ 5, 
6 a 6i, din a a 



(10) 
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FIG. 1: Comparison of difFerent models for the halo mass function originating from Gaussian (dashed) and non-Gaussian (solid) 
initial conditions with the N-body data by PPH08. From left to right, the different panels refers to /nl = 0, 80, 500. Note that 
all non-Gaussian models (with the exception of MR) have been rescaled by the ratio /st//ps, and thus coincide with the ST 
function in the leftmost panel. 



where the new parameter S^, is defined as (5* = 5c\/l — ScS3{a) /3. Since the ST model outperforms the PS one 
in the Gaussian case, it is standard practice to define a new mass function as /mvj /mvj • /st//ps- A further 
modification which has been suggested by [l^ to improve the agreement with numerical simulations is to correct 
the collapse threshold Sc by a factor ^/a = in the expression /mvj//ps- In what follows we will adopt both 

corrections. A similar result was derived by (ssj using a different approach. 

LoVerde et al. (LV) Another way to generalize the PS model is to use the Edgeworth expansion to approximate 
the probability distribution function for the linear density contrast [33| . This gives: 
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LV 
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(11) 



As for the MVJ case, we will use an effective form of this mass function expressing the correction to the ST formula: 
/lv fhv • /st//ps, with the further modification of correcting the collapse threshold Sc by a factor ^/a = in 
the ratio /lv//ps- 

Maggiore-Riotto (MR) Maggiore & Riotto [g^ computed the halo mass function by solving the excursion set 
problem for non-Markovian processes with a path-integral approach, and found 
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Here r(0,x) is the incomplete Gamma function, the additional parameters are a ~ 0.8 and k = an where k 
0.4562 — 0.0040 i?/ with Rf the smoothing scale. Primordial non-Gaussianity affects the function 



where U3 and V3 are given by: 
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The MR function does not need any ad-hoc rescaling. 



12 



d{ai)d{ai) 



{5{al)5{al)5{a^)) 



(13) 

(14) 
(15) 



6 



Lam-Sheth (LS) An extension of the ST model to primordial non-Gaussianity has been recently proposed by 
16411. In this case the mass fmiction is written as: 



/lS ( ^) - /sT 



6 



Ha) 



(16) 



, with f3 = 0.4, 



where H3 (x) = x{x^ — 3), and the mass-dependent coUapse barrier is b{a) = ^/aSc 1 + /? {<J / y/aSc)^'^ 

7 = 0.6, a — 0.7. Note that, in the case of a constant barrier, this mass function reduces to the LV one, if in the latter 
we neglect the term proportional to the derivative of 6*3, which is generally small. 



C. Comparison with N-body simulations 

In Fig. [T] we compare the theoretical mass functions presented above with the N-body data by PPHOS. The halos in 
the simulations were identified using a friends-of-friends algorithm with linking length h = 0.2A, where A is the mean 
interparticle distance. Consistently with PPHOS, for our calculations we assume the WMAP5 ACDM model, with 
parameters h = 0.701, cts = 0.817, = 0.96, f7„ 0.279, = 0.0462, f^A = 0.721. We find that aU the theoretical 
mass functions match the numerical output within ^ 10 — 20% for /nl < 500. Most of the discrepancy originates from 
the fact that the ST mass function underestimates the halo counts from the simulations (left panel in Fig.[T]). Indeed, 
the models are rather accurate in predicting the ratio between the counts in a non-Gaussian model vs. a Gaussian 
one (see also (40l - [43 |). We also consider a fitting formula for the mass function that was computed by PPHOS from 
the very same data plotted in Fig. [TJ Here we rewrite it as 

/PPH 

where we have explicitly introduced a dependence on the threshold collapse density and A, B, C, D are fitting 
parameters which depend on /nl- We use the values from Table 5 in PPHOS. 

In the next section, we will use the mass function to compute the halo bias parameters in the non-Gaussian scenario. 
Given that all the models for n are of the same quality, as a reference, we will only show the results obtained with 
the LV mass function and the PPH fit. 



D + B 



1.686 cr 



exp 



1.6862 0-2 



(17) 



IV. HALO BIAS 

A. Pealt-background split 

We decompose the Gaussian auxiliary potential ip into the (statistically independent) contributions of long- and 
short-wavelength modes: 

<p(q) ^v'Kq) + <p.(q) • (i8) 

Eq. ^ then gives 

= ^1 + hh'p'i - {v^) 

$™ = 2/NL<y3i<^s (19) 

where the dependence on the spatial position is understood. The mixed term $,„ contributes to the short-wavelength 
part but derives from the coupling of Lpi and ips- It vanishes for Gaussian initial conditions. When passing from real 
to Fourier space, the products of two fields become convolutions. Strictly speaking, the terms (piips and tpl would also 
contribute to the long modes due to the mixing of modes caused by the convolution operation. We have checked 
however that these additional contributions are completely subdominant. 

Using the Poisson equation, V^^ = A5 with A — 3ilmHQ / {2c^) , for the density fluctuations we can then write 

61 = SGi{l + 2hi^ipi)+2A-^h-LVipi-Vipi 

Sm = 2hi,{6Gs'Pl+SGlVs)+4.A^'^hi^Vipi-'V(ps 

Ss = Scsil + 2/nl¥'.) + 2A-i/nlV(^, • V^, , (20) 
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where V^i^ = A6g- Notice that: 



2/i 



NL 



1 + 2/nl</?s 



1 + Z/nlV'; 



VnlV(^/ • V(p. 



(21) 



In the spirit of the peak-background spht 0, [sol - fs^ , the short-wavelength modes of the density field collapse to 
form virialized condensations (dark-matter halos) while the long-wavelength ones modulate the halo counts and are 
responsible for large-scale motions. In the non-Gaussian case, however, the halo collapse will also be influenced by 
the long- wavelength modes of ip and Vi^ which contribute to 6m- In a Press-Schechter approach, modulations in 
di will modify the threshold for halo collapse as in the Gaussian case. However, in the presence of non-Gaussian 
fluctuations, the large-scale modes of the pseudo-potential will also alter the statistical properties of the small-scale 
modes in the density field. This provides an additional source of biasing with respect to the Gaussian case. Suppose 
we want to apply the Press-Schechter algorithm to d. The probability that the small-scale fluctuation dg 4- Sm is 
above the collapse threshold 6c — 5i (probability which is obtained by averaging over 5s) would then explicitly depend 
on (5;, Lpi and V(/?i. This implies that the resulting large-scale halo overdensity cannot be proportional to 5i as in 
the Gaussian case (to first order). Rather, in the general case, 5h — bidi -\- fi(pi -\- giVtfi ■ "S/ipi plus higher-order 
terms. The bias coefficients are given by the Taylor expansion of the conditional mass function n[M\6i, tpi, W^pi ■ V</3;). 
Unfortunately, the models for the mass function listed in the previous section have been obtained by averaging over 
the entire Lagrangian volume and have no memory of the cross-talk between large and small scales. We attempted 
the calculation of n{M\5i, (pi, \7(pi ■ V(/Pi) by adopting a Press-Schechter approach and starting from the Gaussian fields 
ip, 'Vip and V^(^ but we could not obtain a closed form due to the complexity of the expressions. 

An approximated model can be obtained assuming that halos form from the highest peaks of 5s* We want to 
implement this requirement in Eq. (j20p . Let us consider what the labels "short" and "long" mean in practical terms. 
The short part of the fields will only include a narrow shell of modes centered around the wavelength corresponding 
to Lagrangian size of the halos. On the other hand, the long part will be formed with all the Fourier modes with 
larger wavelengths. In this case, ips will be closely tracing 5gs — 5. This implies that the high density peaks will 
nearly coincide with the maxima of ips, where V(/?s = 0. In this case. 



5s + 5„ 



1 + 2fni,Lps 



2A-^fniy^i ■ Vipi 

1 + 'ifT^iL'pi 



"^f^-Lfs 



The Lagrangian size of galaxy- and cluster-sized halos ranges between 1 and 10 Mpc. This implies that (5^(y9^)^/^ ^ 
{5f(pl)^^^ and (ipf) ~ i'p'i), since perturbations in the pseudo-potential are nearly scale invariant. Moreover, 
fNhi'pf)^^^ ^ 1 for the values of /nl of physical interest. We thus obtain 



5s [1 + 2/nl<^0 , 



(22) 



i.e. the amplitude of small-scale density fluctuations is enhanced in regions where ipi is large. Therefore, the conditional 
mass function n{M\ipi) can be computed from the unconditional one n{M) by simply multiplying the r.m.s. of the 
density fluctuations by the factor 1 -I- 2/NL'y3i-^ At the same time, we can use the peak-background split to derive 
n{M\di, (fii) by simply replacing 5c with 5c — 5i. 

In each point in Lagrangian space, we can thus define a Lagrangian halo density field as 



Skin) 



n[M, Si{q),ipi{q)] 



1, 



(23) 



where the average can simply be taken as n = n{M, 0, 0). Here it is possible to replace n with / since the proportionality 
factors cancel out, so that we can write more explicitly 



^^(q) 



+ 2/NLipi(q)]'^ 



(24) 



* We only assume that halo formation happens around some of the density peaks. This is different from the approach by |35|. |39|| . in which 

all peaks form halos, and a one-to-one correspondence between them is assumed, 
t Slosar et al. |36| and Afshordi and ToUey |37]] derived a similar expression but notice that ours is written in terms of the non-Gaussian 

density field. 
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We can then expand the perturbations in a Taylor series in terms of both variables 5i and (pi , obtaining 

^.^(q) = EE-S^/(q)^r(q). (25) 

j=0 m=0 

Up to third order in the perturbations, this gives: 

+ ^(6fo'5=' + 3 6^i<5V + 3 6f-2V + &03¥''), (26) 

where all the density perturbations on the r.h.s. are Lagrangian and non-Gaussian. 

Eq. (1^^ implies that not all the coefScients 6^"^ are independent. In particular, all the foj',^ with m 7^ can be 
written in terms of the bj'g. Up to third order we have: 



&01 


— 2/nl 






~ 2/nl 


(-&fo + &20) 




= 4/2^ 




&21 


= 2/nl 




&f2 


= 4/Sl 


(26fo-45,&^o + <5,2 6f„) 


&03 


= 8/nl 





(27) 

It is important to remember that the functional form of the halo mass function accounts for the effect of non- 
Gaussianity on the short wavelength modes 6s, Sm- The bias coefficients bj'Q will then depend implicitly on /nl 
through the shape of the mass function. 



B. Extension to higher-order non-Gaussianity 

If the model for non-Gaussianity is extended to higher order, then Eq. ^ is replaced by Eq. ([1]). If we consider 
cubic corrections, we have the following additional contributions to Eqs. (fT9)) : 

A$™ = 3gNL iff + VI vl) (28) 

which correspond to the following additions in Eqs. ([50]): 

A(5i = 35NL'^?'5Gi +6^"\gNL'^i V<^i • V((5; 

A(5,„ = 35NL^Gi'Ps(2(pi +</3s) + 6gNL^"VsV(p; • V((9; (29) 

A5, = 0, 

where we have already imposed the peak condition. In analogy with the previous section we thus identify the leading 
term as: 

5s + 5yn ~ 5s{^ + ^f-Mh^Pl + 3gNL<y5?) • (30) 

It follows that the r.m.s. of the small-scale density fluctuations will now be altered by a factor (l -|- 2f^i^ipi + ig^i^ipj) 
with respect to the Gaussian case. Therefore, considering ^nl 7^ introduces additional terms in the bias coefficients 
in Eq. (gT]), given by 

A&02 = 65NL^c&fo 

A6f2 = 65NL(-fofo + '5cfo^o), (31) 

while all the other bias coefficients remain unchanged (apart from the modifications due to the implicit dependence 
of the mass function on (7nl, which we do not calculate here). 
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Eq. (pO)) can be finally generalized to an arbitrary order N as 

N 

Ss+S,n^ Ss^jQNLjfi^^- (32) 
i=2 

This equation shows that the leading contribution of each successive order will depend on a higher power of the 
potential, and its effects will therefore be smaller in amplitude. Note that /inl = Qnl4 is the highest-order term that 
can explictly modify the bij parameters (up to third order in the bias expansion), although all the Qnlj will introduce 
implicit dependences by modifying the halo mass function. 



C. Bias from a mass function 



We want now to explicitly calculate the halo bias corresponding to a given mass function. As discussed above, in the 
non-Gaussian case the mass function will also be explicitly dependent on the potential (p, and the halo overdensities 
can now be derived from Eq. (j25p . as a bivariate series expansion in terms of Si and ^pi. Since the effect of the 
short-wavelength modes is taken into account by the functional form of the mass function, we will henceforth drop 
the I indices and use the symbols S and f to denote the long- wavelength parts of the perturbations. 

PS model For reference, we derive the bias coefficients corresponding to the simple PS mass function (see also 





1 


Sc 
^ (72' 


b^o = 


S'c 


3 


b^o = 


S'c 

ct6 


6Sc 



^ (33) 
and show their mass dependence in the right panel of Fig. [21 The remaining coefficients can be obtained using Eq. (j27l) : 

,L f^,26^c 8Sc\ 

- 2(^-10|+4)/2, 

- 2^' M f 

,L _ ./255 20S^ 34Sc 4\ 2 
bi2 - ^1^^-^ + — -^j^NL 

Notice that the "usual" bias coefficients (&io, &20' ^3o) ^^'^ this case independent from /nl. This does not hold in 
general, since an implicit dependence on /nl will be introduced by any non-Gaussian mass function. 

General case To derive the bias coefficients up to the third order, we can repeat the same procedure for any 
other mass function. We have calculated these coefficients for all the mass functions listed in Section Hill finding an 
overall agreement in the trends with mass and with the non-linearity parameter /nl- The analytic form of the bias 
parameters is much more complex than for the PS mass function and we will not write it explicitly. As an example, 
in Fig. we show how the bias coefficients foj'g depend on /nl and halo mass for the LV mass function and the PPH 
fit. 



D. Lagrangian and Eulerian bias 



A model for the formation of the LSS provides a relationship between the density perturbations in Lagrangian 
and Eulerian space. This relation is generally non-local [52] but a local approximation (where x = q) suffices to 
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FIG. 2: Lagrangian bias factors at z = as a function of /nl for a halo mass of 2 • IO^^Mq/Zi (left), and as a function of the 
halo mass, for /nl = 500 (right). The results obtained from the PPH and LV mass functions are shown in both cases. For 
reference, in the right panel we also show the results obtained from the PS mass function, which does not depend on /nl- 



approximately describe the evolution of large-scale perturbations. In this case one writes [51l |5. 

oo 

J=l 

where the Oj's parameterize the evolution of mass-density fluctuations. For the simple case of spherical collapse, we 
have [65|] 

ai = 1; 02 = -17/21; as 341/567. (36) 

Starting from the Lagrangian halo density perturbations Sf^ given in Eq. (|26p we want to use Eqs. and (I55|) to 
write the Eulerian halo overdensity in terms of Eulerian density perturbations. This gives: 

Sh{-x.) = bo + bioS + boiip + 

+ ^{b2o6^ + 2biid(p + bo2^^) + 

+ ^(^^30'53 + 3 621'5V + 3fol2V+^>03<^'), (37) 

where all the density perturbations on the r.h.s. are Eulerian and non-Gaussian (from now on we drop the superscript 
E and all densities will be Eulerian unless explicitely stated otherwise) and the bias coefficients are given by the 
following expressions: 

bio = 1 + 0,1 6fo 

620 = 2(ai + 02) bfo + a? ^20 

^^30 = 6(a2 + 03) b^o + 3 (a? + 20102) b^o + a? ^30 (38) 

bii = 6^i+oi6fi/2 

bo2 = bo2 

621 = (oi +02) 6fi -fo? 6^1/3 
bi2 = 602 + «i ^12/3 

^03 = &03- 
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gio 


= bio 


qii 


= feoi 














920 


= 6io 


921 = 


= 620/2 


922 = fell/2 


923 


= 602/2 








930 


= bio 


931 


= 620 


932 = 611/2 


933 


= 621/6 


934 =612/6 


935 =630/6 


936 =603/6 



TABLE I: Mapping of the bias coefficients in tiie full perturbative expansion. 



Note that Eq. ([57)) differs from Eq. ^ due to the presence of extra terms which are proportional to different powers 
of the Gaussian auxiliary potential tp. Since tp is nearly scale invariant while S has a larger variance on smaller scales, 
the additional terms will only affect the statistics of the halo distribution on the largest scales. Also, (p does not 
evolve with time while S does (according to S{z) = D[z)5{z = 0) at linear order) thus implying that, for a given set 
of bias coefficients, the new terms will become less and less important over time. A third peculiarity of Eq. (j37l) is 
that the halo overdensity can differ from zero also in regions with mean mass density. 



E. Perturbative expansion 



In order to account for the non-linear evolution of mass-density fluctuations in Eq. p7l) we use standard Eulerian 
perturbation theory (see [6^ for a review). We therefore expand the density fields to third order as 5 = 5i + 82 + 5-^ + 
0(54) where 5„ is 0{5i). Each term can be expressed as 

5„(k) = y"^^-^...-|yp-(5£, ^k-^qi^ J„(qi, . . . ,q„)5i(qi) . . .^i(q„) , (39) 

where is the Dirac delta distribution, the tilde denotes Fourier transformation, and the J„ are specific kernel 
functions. On the other hand, since the Gaussian potential ip is the primordial one, there is no need to expand it, 
and it fully coincides with its first-order part Lp = Lpi. 

We can now explicitly rewrite Eq. (1371) up to the third perturbative order as 

(5;i(x) = bo + qioSi + qiiLpi + 

+ q2052 + (3'21^1 + 922<5l'/'l + q2z^\ + 

+ qsoh + q^iSih + q'i2S2Vi + qas^l^pi + qsAipl + q^^Sl + q3,Qip\, (40) 

where, to simplify the notation and facilitate bookkeeping of the terms which will appear in the perturbative expression 
for the power spectrum, we have replaced the biases bij with new coefficients qij . The explicit form of these is given 
in Table m 

Eq. (1301) fully describes the Eulerian halo bias at the third perturbative order in the non-Gaussian case. The leading 
order term, 6h—W + ai^ib(/NL)] ^ + 2/NL&ib(/NL) was already recognized by Dal07, [131 and [36j . 



V. CLUSTERING STATISTICS AND NON-GAUSSIANITY 



Statistical analysis of random fields, such as the Bardeen potential ^'(x), can be performed by studying the irre- 
ducible A''-point correlation functions ($(xi)$(x2)...$(xjv)), or alternatively the A^-spectra 

(2^)35jv(ki, k2, kA,) 6d{\^1 + k2 + ... + )^n) = ($(ki)$(k2)...$(kA,)) . (41) 

For Gaussian fields all odd-order spectra vanish. On the other hand, thanks to Wick's theorem, the reducible even- 
order correlators can be decomposed as products of the power spectrum, 

(27r)3F^(fc)<5z5(ki +k2) = (I>(ki)|.(k2)), (42) 

which, in this case, encodes all the information. If some non-Gaussianity is instead introduced, then higher-order 
statistics become important, such as the bispectrum 

(2^)3s*(ki, k2, ka) 5d{\^i + k2 + ka) = ($(ki)$(k2)$(k3)), (43) 

and the irreducible trispectrum 

(27r)3r$(ki,k2,k3,k4)<5l3(ki +k2 +k3 + k4) = ($(ki )$(k2)$(k3)$(k4)) . (44) 
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Using our simplest model of non-Gaussianity given in Eq. ^ to define the non-Gaussian potential <I> in terms of the 
Gaussian one ip, we obtain: 

P^(k) = P^(fc) + 2/^L J ^P^(g)P^(|k - q|) ^ P^(fc) (45) 

P*(ki,k2,k3) ~ 2/nl [P^(fci)P^(fc2) + (2 cyclic)] (46) 

T<t,(ki, k2, k3, k4) ~ 4/2^ {P^(fci)P^(fc2) [P^(|ki + kal) + P^(|ki + k4|)] + (5 cyclic)} (47) 

where we dropped a sub-leading term proportional to for each of the A'^-spectra (this is why we used the symbol 
of approximate equality). We have checked that the discarded terms are indeed negligibly small. For instance, the 
sub-leading contribution to P* contributes less than 1% of the total for |/nl| ~ 1000 in the fc-range of interest. 

Considering also the third-order term in Eq. ([T]) introduces additional contributions to the power spectrum of the 
Bardeen's potential. The leading-order term can be written as: 



AP$(fc) = Gg^Mk) I ■ (48) 



For Pip{k) (X fc"''^^, the integral above presents an ultraviolet (fc — >■ oo) divergence if rig > 1 and an infrared (fc — !■ 0) 
divergence if < 1 . In general, this is not a problem as the physical process creating the fluctuations will automati- 
cally introduce cutoffs in P^p at small and large wavelengths. For example, cosmic inflation will generate perturbations 
with characteristic sizes comprised between the reheating scale and the present-day horizon However, if AP$ 
is non-negligible with respect to the leading-order contribution (i.e. if 6|(7nl|(<P^) is not much less than unity), 
the results of the perturbative expansion are of limited use unless artificial cutoffs are introduced and the parameters 
of the theory are renormalized. The condition above reduces to |(7nl| lO'' if the currently favored values for the 
amplitude and the spectral index of primordial perturbations are plugged in. Present-day observational limits on 
5nl m, nil therefore suggests that AP$ should contribute at the percent level or less to the power spectrum of the 
potential. Note that in numerical simulations '39], non-physical infrared and ultraviolet cutoffs are introduced by 
the use of a finite volume with periodic boundary conditions. Considering a non-vanishing (7nl also adds another 
leading-order correction to the trispectrum of the Bardeen potential: 

Ar$(ki,k2,k3,k4) ~ 6gNLPAki)PAk2)PUk3) + (3 cyclic) , (49) 

while it does not modify the bispectrum of $ at leading order. 

Linear perturbations in the density at redshift z are related to those in the primordial potential (formally at z — > oo) 
by the Poisson equation 



with 



5i{k) = aik)^{k) (50) 

where the matter growth factor D{z) and the transfer function T(k) have been introduced to account for the linear 
evolution of Si. The function g{z) = (1 + z)D{z) is the linear growth factor for the potential, and g{oo)/g{0) ~ 1.3 
in the currently favored cosmology. Therefore, we can relate the power spectrum of linear density fluctuations to the 
power spectrum of the primordial potential by writing 

Po{k) = Ps,{k) = a'^ik) P^{k) a^{k) P^{k), (52) 



where the last approximation follows from Eq. (|45l) . Similar equations can be written for the three- and four-point 
correlators of the linear density perturbations, which we will label Bq and Tg respectively, by combining Eq. (|50p with 
Eq. dUl) and gll). 



In multi-field inflationary models the non-Gaussian contribution to the trispectrum may scale independently from the bispectrum. For 
this reason, the factor /^j^ in Eq. I|47|l is sometimes re-labeled tnl, and treated as an independent parameter. Observational constraints 
on the trispectrum of the Bardeen's potential should then discriminate between such models and the simplest inflationary scenarios. 
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The fact that the Bardeen's potential decays with time proportionally to g{z) implies that the actual values of the 
coefficients Qnlj depend on the cosmic epoch at which Eq. ([T]) is applied (see Section 2.2 in [4l|). Here we apply it at 
early times (which is sometimes called the "CMB convention" ) while other authors use the fields linearly extrapolated 
at z = (the "LSS convention"). In general, 



nLSS _ ^CMB 



5(0) 



(53) 



so that — 1-3 /nl'^ ^^^d — I-'^Snl'^- This conversion factors should be taken into account when comparing 
papers using different conventions. 

VI. POWER SPECTRA 

We are now ready to calculate the two-point statistics of the LSS arising from non-Gaussian initial conditions. We 
compute the halo- halo power spectrum and the halo-matter cross spectrum as follows. First, we take the Fourier 
transform of Eq. (|40|) and build the corresponding two-point correlators {Sh(ki)Sh(k2)) and {Sh(ki)S{'k2)) . These are 
composed of many pieces and we only consider terms up to the fourth perturbative order. For instance the leading 
contribution to halo-halo spectrum (second order in terms of the perturbations) is composed of 3 terms, the third order 
correction is made of 8 pieces (of which one identically vanishes because is a Gaussian field) , and the fourth-order 
one contains 30 terms (14 of which are obtained multiplying a linear perturbation by a third-order one - indicated 
by the subscript (13) hereafter - and 16 are originated by the product of two second-order terms - subscript (22) 
hereafter) . To proceed we then: (a) use Eq. ((39|) and write the density perturbations of order n > 1 as convolutions of 
n linear perturbations and a kernel J„ ; (b) express the linear density perturbations in terms of the potential using 
the Poisson Eq. (|50|) : (c) take the ensemble averages by using the expressions for the power spectrum, bispectrum 
and trispectrum given in Eqs. PSI) . (IITI) . P5|) . and ([5^ . While ((^(ki)(^(k2)<^(k3)) — 0, attention must be payed to 
the mixed terms in $ and ip as: 

(I>(ki)l>(k2)l>(k3)) ~ ((^(ki)l>(k2)l>(k3)) + (2 eye.) ~ ((^(ki)^(k2)l>(k3)) + (2 cyc.) cx /nl (54) 

where equalities only hold at leading order in ip" (i.e. ip'^) as the central correlator also contains a sub-leading term 
proportional to /nl^nl (which scales as ip^) and the leftmost one some terms proportional to /^li /nl^nl (both 
scaling as </?^), and /nl^nl V^)- Similarly, to leading order in cp" (i.e. (p^), 

(I>(ki)l>(k2)l>(k3)l>(k4)) ~ fLTA+g^LTB (55) 

(^(ki)l>(k2)l>(k3)$(k4)) + (3 cyc.) ~ 2 J^^Ta + 3 gNLTe (56) 
((^(ki)^(k2)l>(k3)$(k4)) + (5 cyc.) ~ f^^TA + 3g^LTB (57) 
((^(ki)(^(k2)(^(k3)$(k4)) + (3 cyc.) ~ ^nlTb (58) 

where TA(ki, k2, k3, k4) and rB(ki, k2, k3, k4) are defined in Eqs. (|T7)) and 

A. Results and comparison with N-body simulations 

Matter If we set bio = 1 and all the other bias coefficients to zero, we obtain an expression for the power spectrum 
of mass-density perturbations which coincides with the result by P"^"^{k,z) = D'^{z)Pii{k) + D^{z)P"2"'{k) + 
D^{z) [P2™"(fc) + P^'^ik)], where 

P™'"(fc) = Po(fc) 

Pink) = 2 j ^4^)(q,k-q)Po(-k,q,k-q) 
P2r{k) = 2 j ^[j(^)(q,k-q)]'Po(g)Po(|k-q|) + 

^2^)6 -^2 (P,k-p) J;^ '(q,-k-q)To(p,k-p,q,-k-q) 
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FIG. 3: Left: Deviation of the matter power spectrum in models with different /nl (and giNL ~ 0) from the Gaussian case at 
z = (top) and z = 1 (bottom). The lines indicate our one-loop calculation for different values of /nl while points with error 
bars correspond to the N-body simulations by PPH08. Right: Halo- matter cross spectrum at 2: = for a narrow bin of halo 
masses centered around M = 2 • 10^'^ Mq / h (top) and at z — 1 for M — 5 ■ IO^'^A/q/Zi (bottom). The solid and dashed lines 
have been obtained using our model with the bias parameters from the LV and PPH mass functions, respectively. The dotted 
lines indicate the model by Dal07, while points with error bars correspond to the simulations by PPH08. 
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\k) = 6 1 |!3_j(^)(k,q,_q)Fo(5)F„(fc) + 



(27r)3 

+2 J ^^|^ji'^(p,q,k-p-q)To(-k,p,q,k-p-q), (59) 

(s) . . 

and Jn indicates a kernel which has been symmetrized with respect to its arguments. We have checked that the 
two-loop contributions proportional to Tq are, in general, negligibly small and will not be considered hereafter if 
ffNL = 0. In the left panel of Fig. [3] we plot the ratio between the matter power spectra originating from a non- 
Gaussian model (with /nl 7^ and gNL = 0) and from Gaussian initial conditions. We consider several values of /nl 
and we compare our analytical results with data from the N-body simulations by PPHOS at both redshift and 1. 
Primordial non-Gaussianity alters the matter power spectrum at the few percent level for k < 0.2 h Mpc~^ and these 
deviations arc remarkably well reproduced by the one-loop corrections. 

Halos The halo-halo power spectrum (and similarly the halo-matter cross spectrum) deriving from our multivari- 
ate biasing scheme can be written as 

P'^^ik, z) = D\z) P^iik) + D\z) P^iik) -f D\z) [pU{k) -f P^ik)] , (60) 

where the superscript j indicates either matter (m) or halo {h) fluctuations. The full expressions of the different 
terms are lengthy and we report them only in the Appendix. We highlight that our calculation reduces to: (a) the 
linear result by Dal07 if we only consider the leading-order terms (and further assume that the mass function does 
not depend on /nl); (b) the usual one-loop Gaussian expression derived e.g. by ^57|] if we set /nl = 0, and (c) the 
non-Gaussian result by [43| if we ignore the terms which are proportional to the potential perturbations f in our 
multivariate biasing scheme. Notice that for fc — we obtain P^^{k) oc PQ{k)/a'^{k) and P^™'{k) ex. Po(fc)/a(/c). 
In the right panel of Fig. [31 we test our theoretical predictions for the halo-matter cross spectrum as a function of 
/nl and for ^nl = (solid lines) against the N-body data by PPHOS. The bias factors in the models have been 
calculated from the LV and PPH mass functions. We consider halos with mass M ~ 2 • 10^'^ Mq/H at z = and 
M ~ 5 • lO^^MQ/h at z = 1. We have chosen two different mass bins to keep the number of halos at each redshift 
large enough to avoid substantial shot noise contamination in the simulations. Our analytical results are in very good 
agreement with the simulation data for the whole range of /nl and up to scales k ^ 0.2 h Mpc"^. Note that our 
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spectra differ from the linear result by Dal07 (over-plotted with dotted lines) both on large and small scales. The 
small-scale departure is due to the non-linear growth of perturbations which we take into account up to the third 
perturbative order. The large-scale discrepancy is discussed in detail in the next subsection. In Fig. |4]we show how 
the cross spectrum depends on the halo mass at three different wavenumbers and for two redshifts. Independently 
of halo mass, at z = 1 our model accurately matches the outcome of the simulations for k < 0.2 h Mpc^^. Likewise, 
at z = 0, the theory agrees well with the numerical data on the largest scales while it tends to over-predict the cross 
power for fc > 0.1 /i Mpc~i and M < IO^Mq/Zi. 

We have checked that, for ^nl — 0, the terms proportional to the trispectra of the potentials $ and ip are generally 
subdominant even when they generate contributions to the halo power spectrum which diverge as fc — 0. For instance, 
the trispectrum contribution to the term (SfSf) in the halo-halo power spectrum scales as Po{k) /a'^{k) as like as the 
leading term. However, for |/nl| < 500, this correction contributes at most at percent level and only on very large 
scales. Also the trispectrum contribution in {Sf Si) which scales as Po{k) / a{k) is subdominant (^ 1%) for both the 
halo-halo and the halo-matter cases. Similar conclusions can be drawn for the terms including ip: we have checked 
that the contributions arising from averages where one or more are replaced by ip are subdominant. 

The effect of gNL The situation becomes more complicated if we consider also the third-order term in Eq. ^ 
with realistic values of (^nl- In this case, there will be new contributions to the halo and halo- matter power spectra 
coming from two sources: the trispectrum gets the additional term of Eq. ()49|) . and the bias factors become altered 
as described in Eq. pil) : each of these modifications is linear in ^nl- 

Considering first the effect of AT$, we have found that this adds a (negligible) constant contribution to (SfSf) but 

generates another term which scales as Po{k) / a{k) in (SiSi). The latter can become the dominant contribution on 
large scales and for high values of ^nl ^ 10^. In the limit fc — > 0, this two- loop term (whose full expression is given 
in Eqs. (|82l88p in the Appendix for the halo and halo-matter cases) reduces to 

a(k) 



1 . _4,m^ .m^o(fc) 



where we define 



PiP (k) ^ :.<?NL&30^'(i?)S3(i?)^(X.gNLfc"-^ (61) 

2 a(K) 

^3(i^) - / ^ ^ ^ «(|p + q|) Wm WipR) , (62) 

such that 5*3 — f^j^ E3. In our formalism, P'^^'^\k) corresponds to the leading contribution found by [39| . who used it 
to derive observational constraints on ^nl (assuming /nl = 0)- In the limit of high peaks, the scaling 610^30 ~^ ('^c/c)'^ 
is recovered. Note, however, that the amplitude of this term depends on the adopted smoothing scale R first introduced 
in Section ini and further discussed in Section IVl CI 

Second, we look at the effect of A&. As shown in Eq. ([3T|). only the coefficients &02 and &12 are altered by gNL- 
In the halo- matter case these new terms are subdominant with respect to the trispectrum contribution of Eq. (|6ip . 

However, in the halo-halo spectrum, the leading term for fc — > is f(23)(23) ^ °c which scales as 

^(2'3)(23) -> \ bl, al{R) ^ (X 4l , (63) 



where cr^ = / dqq^ [Pi^{q) / [q)] W''^ {qR) / {2ti^) . Because of the quadratic dependence on gNL, this term dominates 
on very large scales for high values of ^nl and small /nl- Its dependence on the smoothing radius R is very weak 
because the potential </? is nearly scale invariant. Note that the terms in /nl5nl are generally subdominant. 

B. Bias and asymptotic behavior on large scales 

Let us define the effective bias function 

P^"'(fc,/NL) 

^cff(fc,/NL)^ ^_(^^^^^) , (64) 

and compare it with the standard Gaussian bias by introducing the bias deviation 

A&(fc, /nl) = h,fi{k, /nl) - 6cff(fc, 0). (65) 
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t„, = 500, z = 1 




FIG. 4: The halo-matter cross spectrum as a function of halo mass at z = (left) and 2 = 1 (right) is plotted for three 
different values of the comoving wavenumber k (in h Mpc"'^). Colored solid and dashed lines indicate the results of our 
perturbative calculation using the LV and PPH mass functions, respectively. Data points with error bars correspond to the N- 
body simulations by PPH08. The thin lines show the linear theory by Dal07 (black), also corrected with the factor /3 (magenta) 
introduced by PPH08 (see Section TVI Bl for further details). At 2: = structure has evolved further into the non- linear regime, 
so that the range of validity of both the linear and one loop theories is reduced. 



In this section we will only consider the case ^nl = 0. In the limit fc — and for large R {a\ ^1), the dominant 
contribution to the halo power spectrum is given by the tree-level term, and we thus obtain 

A6iincar(fc) = 6io(/nl) - &io(/nl = 0) + 2/nl5c [&io(/nl) - I] I a{k) . (66) 

The scale-dependent non-Gaussian correction is proportional to the factor &io — 1 as originally shown by Dal07, 
although there is also an additional scale-independent correction, due to the fact that in our model 610 is a function 
of /nl (similar conclusions have been reached by [s^, [38l - l4l| following different approaches) . In the simplest model 
by Dal07 the scale-independent term is missing: 

A&Dai07(fc) = [feio(/NL) - l]/a(fc) . (67) 

In the left panel of Fig. [SI we plot Ab{k) for /nl — 500 and test the different models against the N-body simulations by 
PPH08. We can see that considering only the scale-dependent term as in Dal07 does not match the simulations very 
well, since, contrary to the N-body data, A6 cannot change sign with increasing k (see also |41]). The agreement vastly 
improves if we use Eq. (|66p with the bias parameters computed from a non-Gaussian mass function (LV) as this adds 
a constant negative shift (see the left panel of Fig. [2]) to the bias deviation. Considering the full calculation to third 
perturbative order further improves the agreement with the simulations for k > 0.1 h Mpc~^ up to a maximum value 
of the wavenumber which depends on the adopted smoothing scale for the perturbative calculations (see Section [VI CI 
for further details). To highlight the importance of the non-linear and scale- independent corrections, in the right 
panel of Fig. [5] we plot the ratio 

A6(fc,/NL) . . 

A6Dal07(A;,/NL) ^ ^ 

for several values of /nl and using both the LV and PPH mass functions. We consider halos with mass M = 
2 ■ lO"/i-iM0 at z = 0. The simulation data by PPH08 are in good agreement with our third-order calculation, 
while the linear (Dal07) model cannot reproduce them on scales k > 0.02 h Mpc^^. The impact of primordial non- 
Gaussianity on the halo bias is further explored in the left panel of Fig. [51 where we show the bias deviation as a 
function of /nl at three selected scales, again for halos with M = 2 ■ lO^'*/i~^M0 at z = 0. Note that, contrary to 



17 




_U I I I I 

0.01 0.1 
k [h Mpc-'] 




k [h Mpc-'] 



FIG. 5: Left: Change in the effective bias A& for fixed /nl = 500, M = 2-10^^ Mq/H, at ^ = 0. Different models are compared 
with N-body simulations: the simple AbDaioT (dotted), Abuncar (dashed) using the LV mass function, and the full one- loop 
theory (solid), which yields the best match. Right: Fractional difference between the one-loop prediction for Ab (with LV and 
PPH mass functions) and the Dal07 linear theory, compared with the simulations, for different values of /nl, at the same mass 
and redshift. 




FIG. 6; Left: Bias deviation at 2 = as a function of /nl at three different wavenumbers, and for M = 2 ■ 10^* Mq/H. The 
results for LV (solid) and PPH (dashed, within its range of validity) mass functions are shown. We also plot the prediction for 
the linear (Dal07) theory. Right: As in the left panel but as a function of fef and for /nl = 500. The normalization factor 
a/r = a/{2A) is chosen to reproduce figure 11 in PPH08. 
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what predicted by the hnear (Dal07) model, the relationship between A6 and /nl is non-linear, as first observed by 
PPH08 in their N-body simulations, and this is now fully explained by our perturbative calculation to third order. We 
finally study the mass dependence of the bias deviation by showing, in the right panel of Fig. |6l how Ab changes as a 
function of the first Lagrangian bias coefficient fofg. We consider three values of the wavenumber in the quasi-linear 
and mildly non-linear regime, /nl = 500, and z = 0. In the Dal07 model, Ab depends linearly on and the relation 
between these two quantities is independent of k. This is indicated by the solid black line which departs more and 
more from the simulation data with increasing k. In order to better describe the N-body results, PPH08 presented a 
fitting function, /nl), in the form of a multiplicative (scale-dependent) correction to the Dal07 model (magenta 
lines). Note that our third-order calculations match well the numerical outcome. Some discrepancy is noticeable for 
k = 0.2 h Mpc~^ and bfg < 1, where non-linear effects become more important (due to the small first bias coefficient) 
and perturbation theory becomes less accurate. The difference between the bias deviations obtained with the LV and 
PPH mass functions at large halo masses emphasize the need for accurate parameterizations of the halo counts for 
the rarest objects. 

Our result in Eq. (|66p differs from the perturbative calculations based on the univariate local bias by [i^, where 
the scale-dependent part of the bias deviation was found to scale as 620 times the variance of the mass density field. 
Strictly speaking, for /nl 7^ and k our result gives P'"^{k, fyh) {boi + 262o/NLO'fl)Po(A:)/Q!(fc) but the term 
proportional to 620 is suppressed by smoothing on the scale R (which is necessary to truncate the bias expansion 
at third order in a meaningful way). Note that, using the PS expression for the bias parameters in the limits of 
high peaks, Sc/a > 1, this reduces to P^'^{k, f^i^) — >• 2/nl(<5c/ct^)[1 -I- {a^/ ij'^)]Po{k) / a{k) . In this case, the two 
contributions are identical if R is chosen to be the Lagrangian radius of the halos (i.e. R = Rf) and the same 
window function is used to compute a and in the calculation of the perturbative power spectra. In general, however, 
using the Lagrangian radius of the halos gives a of order unity and this is too large to allow the truncation of the 
bias expansion at a finite order. For this reason in our calculations we use an < a and the smoothing-dependent 
contribution proportional to 620 is subdominant.'* 

As highlighted in Eq. p7|) the leading order for (5/i in our multivariate biasing scheme includes a term proportional 
to and this generates the scale-dependent correction in Ab. The proportionality with 620 found by other authors 
derives from the assumption that a local deterministic bias scheme holds true also in the presence of non-Gaussian 
perturbations. In this case, for k — > 0, second-order terms dominate over the tree-level contribution to the power 
spectrum which casts some doubts on the validity of the perturbative expansion. In the left panel of Fig. [7] we show 
the difference between our multivariate approach and the standard local bias. The asymptotic scale dependence 
P'*™(fc) oc or^ik) Po{k) cx for /c — is recovered in both cases, but the amplitude of the diverging term in the 

standard local bias model depends on the smoothing length that has to be introduced to cure the ultraviolet divergence 
of the mass variance. In Fig. [7| we use a smoothing scale of 10 Mpc for both models and the asymptotic term 
deriving from the standard local bias is strongly subdominant with respect to the correction given in Eq. (j66p . As 
discussed by [i^, the result of the univariate local model depends strongly on the smoothing scale, as it is proportional 
to (T^(i?). For instance, using R = 2h^^Mpc boosts the amplitude of the scale-dependent bias (see Fig. [7] (left)). This 
is why only the multivariate model can reproduce the results from N-body simulations by PPH08 without tuning 
additional parameters. The main practical advantage in this case is that the bias parameters can be predicted from 
a model for the mass function, while the results from the standard local bias can only be used after "renornializing" 
the bias coefficients [66| and using them as fitting functions. However, even though one can play with the parameters 
of the theory to fit some data, one should not forget that the physical origin of the scale-dependent bias is that 
large-scale fluctuations in Sh trace perturbations in ip and this is not accounted for by the standard local bias model. 
The differences between the models are further highlighted in the right panel of Fig. [71 where we plot the ratio of 
the halo-matter cross spectra obtained with different approximations with respect to our full non-Gaussian one-loop 
calculation, at /nl = 500 and at z = and 1. This figure summarizes all the conclusion we have reached in this 
section: (1) the univariate local biasing assumption yields the correct fc-dependence but the wrong amplitude of the 
spectrum on large scales; (2) the linear approximation in Eq. (1661) lacks small-scale power; (3) the simpler model by 
Dal07 also features a scale-independent offset in the effective bias. Similar results can be obtained for the halo-halo 
power spectrum, for which the asymptotic scale dependence is P'^^{k) oc a~'^{k) Po(fc) « fc"''"'* for fc — >■ . 



S It is interesting to see what alternative approaches find regarding this discrepancy. For instance, Matarrese & Verde [35j computed 
the two-point correlation function of regions where the density exceeds a high threshold 5c S> O" and found that, for large separations, 
the non-Gaussian correction scales as ((5^/(7^) ^3(xi , xi, X2) with ^3 the three-point correlation function of the mass density. In Fourier 
space this coincides with the high-peak limit of our result but, provided that = a^, it also matches the result by [4^ . This happens 
because for high peaks both (5c(6fQ)^ and bfgfelo are proportional to <5^. The same ambiguity applies to the higher-order calculations in 
Desjacques & Seljak [39j . 
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FIG. 7: Left: The standard univariate local bias prescription (local) is compared with our multivariate scheme for the halo 
bias (LV full). We consider the same halo masses as in Fig. [3] (right) and we plot the halo- matter cross spectra deriving from 
the two bias models at 2 = 0, 1 and with a smoothing scale R = 10h~^ Mpc. In the Gaussian case (/nl ~ 0) the cross spectra 
coincide, but they are very different for /nl = 500. While the asymptotic scale dependence for A: — >■ is recovered in both 
cases, only the multivariate model can reproduce the results from the N-body simulations. We also show the dependence on R 
in the univariate model, by overplotting the result with R = 2h~^ Mpc. Right: Ratio between the halo-matter cross spectra 
obtained with different approximations and our full perturbative calculation, using /nl = 500 and for the same redshifts and 
halo masses as in FigO Assuming univariate local biasing (red, solid) severely under-predicts the large-scale power. On the 
other hand, the linear approximations in Eqs. (|66|) (green, dashed) and (|67p (blue, dotted) lack small-scale power. Notice that 
the Dal07 model also features a constant offset on large scales due to the missing scale-independent correction discussed in the 
main text. 



C. On the smoothing and other perturbative approaches 

The series expansion in Eq. ([2]) only applies when Sh and S have been smoothed on a scale R. The reason is twofold. 
First, the locality of the bias is expected to degrade progressively when smaller and smaller scales are considered. 
Second, we truncate both the bias and the perturbative expansions to finite order which is a good approximation only 
if the neglected terms give a small contribution. This requires that typically i5 ^ 1, i.e. R ^ 5h~^ Mpc. 

As explicitely indicated in the Appendix, in this paper we have used a Gaussian kernel W{kR) = exp[— (fci?)^/2] 
with R = 10 Mpc to smooth the evolved fields 5i, 62 and ^3 before applying Eq. ([2]). An obvious consequence 
of this procedure is that also the resulting halo power spectrum is suppressed on scales k > R~^. A common 
prescription to lessen this effect and extend the theory to (slightly) larger wavenumbers is to divide out W^{kR) from 
the perturbative result for P'^^{k) and P'""(fc), as introduced by [67|]. Here we have followed this approach to consider 
wavenumbers up to /c ~ 0.3 /i Mpc^^. 

Some of the one-loop corrections to the halo power spectrum present ultraviolet divergences that are automatically 
cured by using a finite value of R. However, some of these integrals give rise to scale-independent contributions for 
A: — > whose amplitude depends on R, as shown already by [STj . This is somewhat unsatisfactory, since it makes the 
results dependent upon a non-fundamental quantity, and it is amongst the reasons which have led to the application 
of renormalization techniques to the theory, often borrowed from other areas of physics. The existin g ap proaches, 
as recently reviewed by [68|, include the renormalized perturbation theory 69, 70], the closure theory [7l[, the time 
renormalization group flow model [t^ , and the renormalization group perturbation theory [66l [t^ . Most of these 
approaches do not include a bias model and just apply to the matter density field. The renormalization of the bias 
parameters, included in some of the models, makes the theory free from any undesired dependence on the smoothing 
scale. This can be achieved by grouping different perturbative terms together and relabeling some parameters to 
include the smoothing-dependent factors, a procedure which is not uniquely defined. An altogether different approach 
which does not need such an operation is the Lagrangian resummation theory by (sgI . [tJ . 
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On the other hand, SPT has the advantage of remaining a fully predictive theory, where the bias coefficients can be 
calculated as a function of halo mass. Furthermore, the choice of the smoothing scale R is not completely arbitrary, 
but confined to a rather narrow range around i? ~ 10 Mpc. Indeed, the smoothing needs to be i? ^ 8 h~-^ Mpc 
in order not to break the validity of the perturbative expansion in a significant fraction of the volume (a <^ 1 for 
matter and a <C 620/^10 for halos) and, on the other hand, R needs to be as close as possible to this limit if we want 
to prevent the smoothing from wiping out the non-linear corrections at the wavenumbers of interest. 

We have checked numerically that different choices of the smoothing scale R within a reasonable range larger than 
the Lagrangian size of the halos do not affect our results significantly. Notice that, for non-Gaussian perturbations, 
the A:-independent - but i?-dependent - terms arising in SPT on large scales are less important, since the halo power 
spectrum grows with decreasing k. 



VII. BISPECTRA 



The leading contributions to the halo bispectrum from non-Gaussian initial conditions of the local type have been 
recently computed in the framework of the local bias model given in Eq. ([2]) [43 - l46j : 



B^(ki, k2, kg) = blBs{ki,k2, kg) + 



Psiki)Psik2) + (2 cyc.) + - 



q 



r5(q,ki-q,k2,k3) + (2 eye.) 



(69) 



2 J (27r)3" 

Using Eulerian perturbation theory to follow the growth of density perturbations gives up to fourth order 

B5(ki,k2,k3) ~ Bo(ki,k2,k3) +2^^2(ki,k2)Po(fci)/'o(fc2) + (2 cyc.) , (70) 
where the second term is generated by non-linear gravity while 

Poih)Poik2) 



Bo(kl,k2,k3) ~ 2/] 



NL 



"(fca) ■ 



a{ki) a{k2) 



2 cyc. 



(71) 



is the linear matter bispectrum due to primordial non-Gaussianity. Similarly, the term between square brackets in 
Eq. (p^ reduces to 



Po(fci)Po(fc2) + (2 cyc.) + i 



d3 



q 



To(q,ki-q,k2,k3) + (2 cyc.) 



(72) 



In full analogy with the power spectrum calculation discussed above, it is straightforward to show that our new 
expansion of 5h in terms of both 5 and gives rise to many additional contributions. The most compact form is 
obtained by writing the expansion of the product of three 5h evaluated in real space at three different locations. 
For the term which generates the contribution to the halo bispectrum which is proportional to the linear matter 
bispectrum, we have: 



bi5i5i5i feio(5i(5i(5i -I- &io^oi<^'5i(5i + (2 cyc.) + binhlitpipS + (2 cyc.) + hl^tpiptp . 
On the other hand, for the term accounting for the non-linear evolution of the density, one finds: 

hl5i6i52 + (2 cyc.) ^ blQ5i5i52 + (2 cyc.) + bl^hQiip5i52 + (5 cyc.) + 610601 w52 + (2 cyc.) 
Finally, for the source of the term between square brackets in in Eq. (1691) . we get: 



(73) 
(74) 



h)lb25i5i5l + (2 cyc.) ^ ]^b%b2Q5i5i5l + (2 cyc.) + bl^hii5i5i{ip5i) + (2 cyc.) -I- h]wh2oboi5iv5l + (5 cyc.) -I- 

^&?0^02^i'5i<y5^ + (2 cyc.) -I- 6io6oi6ii^i^(.^(5i) -I- (5 cyc.) + h3l^b2oVvSl + (2 cyc.) + 
^&io&oi&02'5i<y5</5^ + (5 cyc.) + bl^biiLpLp{Lp5i) + (2 cyc.) + h3l^bQ2VVv'^ + (2 cyc.) . (75) 



The halo bispectrum can be computed by Fourier transforming the expressions above. For instance, from Eq. ([73 
we obtain: 



B,,(ki,k2,k3) 







h , 




bw ^ 


^01 


610+ . 




a[k2) 






a{k3)_ 



Bo(ki,k2,k3) , 



(76) 
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FIG. 8: Theoretical bispectra at tree level in the equilateral configuration: B'^"'"'(k, k, k), B^""^{k, k, k), B'''"^{k, k, k), 
B^^^{k, k, k) for two values of /nl- In the Gaussian case (/nl = 0) the tree-level bispectrum vanishes. 



while assuming a local-bias scheme would have given i?/i(ki, k2, ka) ~ i3o(ki, k2, ka). The mixed matter-halo 
bispectra can be obtained in a similar way, and their expressions differ from Eq. (|76p only for the presence of a 
reduced number of scale-dependent bias factors. We plot the equilateral configuration of the bispectra in Fig. |8l for 
two values of /nl- Since we are only considering the tree-level contribution, the Gaussian bispectrum is vanishing. 
Note that our leading-order result of Eq. ([75)) can be reproduced by taking the analogous formula obtained with the 
univariate local bias and simply replacing &io with the scale-dependent bias 6io -I- h^i /a{k). More complex equations 
relate the scale-dependent terms obtained by taking the Fourier transform of Eqs. (l74t and (|75|) . We will not discuss 
them in detail here. 

It is important to notice that the scale-dependent bias changes the shape dependence of the halo bispectrum. 
The tree-level term Bq is dominated by the squeezed configurations (where one of the wavevectors is small) and 
the /NL-dependent term hoi /a makes it even more so. Comparing our results with the figures in [i^ suggests that, 
in strict analogy with the result for the power spectrum, the terms proportional to boi/a give by far the dominant 
contribution to the bispectrum on large scales. The shapes of the bispectrum parts coming from the non-linear growth 
of perturbations and second-order biasing are very different [45l . l46j and this makes the bispectrum a promising tool 
to measure /nl- 

We can see that the expanded form of the bispectrum deriving from Eqs. (|74p and ([75]) depends on all the non- 
Gaussianity parameters in a non-trivial way: all terms involving the matter bispectrum Bq, i.e. all terms involving an 
average over three (5's or (/s's bring in a linear dependence on /nl- Then, all terms involving the matter trispectrum Tq, 
i.e. averages over four 5's or ip^s, have two contributions, depending on ^nl and /^l (o'" '^^ general tnl) respectively. In 
addition to this, we have additional dependences on /nl and (;nl implicit in the bias factors, as described in Section 
IIVI This complex shape and scale dependence of the bispectrum offers an unique opportunity to simultaneously 
constrain /nl, ffNL and tnl and thus distinguish between infiationary models. We will investigate this in the near 
future. 



VIII. CONCLUSIONS 



We have studied the growth of structure from non-Gaussian initial conditions of the local type. In particular, 
we have shown that the spatial distribution of dark-matter halos is naturally described by a multivariate local bias 
scheme where the halo number density depends on the underlying values of the density field S, the auxiliary Gaussian 
potential ip, and (possibly) also on its gradient. This bivariate local approach can be equally interpreted as a non-local 
description in terms of 6 only, since ip and S are related by the Poisson equation. Adopting the peak-background split. 
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some common parameterizations of the halo mass function, and a local model for the evolution of large-scale density 
perturbations, we have derived the coefficients of this multivariate expansion as a function of the halo mass and of 
the parameters quantifying the level of primordial non-Gaussianity. 

Using SPT to approximate the non-linear growth of density perturbations, we have computed the halo power 
spectrum and the halo-matter cross spectrum up to the third non-vanishing perturbative order. For unbiased tracers 
our result coincides with the matter power spectrum presented by [1^. However, in the most general (biased) case, 
it differs from what is obtained adopting the standard local bias expansion in terms of the density field [43l446l | . The 
most remarkable feature is that the scale-dependent bias first discussed in Dal07 appears at leading order in our model 
for the power spectrum. This is because in our multivariate biasing scheme halo fluctuations on large scales trace the 
Gaussian potential Lp rather than 5. However, our model reduces to the usual univariate case on larger scales, where 
the variance of density fluctuations is much larger than that of the potential. Note that both the multivariate and the 
univariate models predict that the dominant contribution to the halo power spectrum scales as fn'L Po{k) / a{k) for 
A: — >■ 0. However, for the standard univariate biasing, the amplitude of this term is given by a badly behaved integral 
which strongly depends on the assumed smoothing scale. Renormalization of the second bias coefficient (which should 
then be treated as a fitting parameter when comparing the theory to observation or simulations) is unavoidable in 
this case, while it is not needed in our model. 

We have then tested our results against the N-body simulations by PPH08, finding excellent agreement for both 
the matter and the halo two-point functions. Focusing on the scale-dependent bias generated by primordial non- 
Gaussianity, we have shown that our model accounts for the discrepancies previously found between the predictions 
by Dal07 and the outcome of numerical simulations j40l - |43| . Corrections to the simpler model by Dal07 arise for 
two main reasons: (a) the bias coefficient 5io depends on fm, due to the fact that the shape of the mass function 
is altered by primordial non-Gaussianity (see also [H, HI, liy), and this adds a scale- independent offset to the bias 
deviation A6; (b) considering perturbation theory up to third order generates numerous additional corrective factors 
that become important on intermediate and small scales. With our one-loop calculation of the power spectrum, the 
range of validity of the theory extends up to scales k ^ 0.1 — 0.3 h Mpc~^ depending on halo mass and redshift. 

We have also shown how our calculations can be extended to include higher-order terms of primordial non- 
Gaussianity, for instance by considering a non- vanishing primordial trispectrum proportional to the parameter ^nl- 
In this case, the halo power spectrum includes an additional contribution proportional to g^i^Po^k) / {k) which, 
depending on the values of /nl and ^nl, may be dominant on the largest scales. This term only appears in our 
multivariate expansion and originates from the bias parameter 602 which includes a correction proportional to (7nl- 
On the other hand, in agreement with [s^, we have found that both the halo-halo and halo- matter spectra acquire 
a dependence on (7nl from the trispectrum of the linear density field. This term scales as g-Ni^Poik) / a{k) but its 
normalization depends on the assumed smoothing scale and cannot be robustly predicted by the theory. 

Finally, we have calculated the halo bispectrum deriving from our multivariate biasing scheme. At tree level, our 
result corresponds to the usual bispectrum deriving from Gaussian initial conditions but where the linear bias 610 is 
replaced by 610 -I- hf)i/oi{k). This is different from what has been found assuming univariate local biasing [i^ l46j . 
Therefore the analysis of three-point statistics represents a promising tool to test the different biasing schemes against 
observations. Also, the complex shape and scale dependence of the halo bispectrum offers an unique opportunity to 
simultaneously constrain /nl, ffNL and tnl and thus put entire classes of inflationary models under scrutiny. We will 
explore this in more detail in a forthcoming paper. 
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Appendix: Complete analytic expression of the power spectra 

We list below the non-vanishing contributions to the halo-halo and halo-matter power spectra up to one-loop in 
perturbation theory. These have been obtained by smoothing the evolved density perturbations with the filter W{kR), 
so that S{k) — >■ d{k, R) = S{k) W{kR). We highlight with the label "Local" the terms that are also present if we use 
the univariate local bias approach as in [isj . Similarly, the contributions that do not vanish in the Gaussian case are 
marked with the label "Gauss". We have dropped the two- loop terms which arise from the trispectrum To, as they 
generally give negligible contributions, with the exception of the term (k, R), which is important in the case of 
large ijnl and small /nl- Some of the integrals below present an infrared divergence if ~ 1, like for instance the 
term f(23)(23)- ^^^^ case, we introduce a cutoff in Po{k) for k < kn = ^/Rh where Rh = c/Ho [75]. 

Note that only the halo- matter cross spectrum has been compared to the N-body results by PPH08, as in the 
simulations the halo-halo spectrum is more strongly affected by shot noise. 



1. The halo-halo spectrum 

The full halo-halo power spectrum at one loop is 

P'^\k, z, R) = D^{z) Plfik, R) + D^{z) Plfik, R) + DHz) [Plf{k, R) + P^^iK R)] , 

where: 

Pil^{k, R) is the sum of the following terms: 



Local, Gauss Pll^o-^f^j^g-^ik, R) = bi„ Po{k)W^kR) 

(10)( 



^P(\o)(ii){k,R) = 2feio6oi^^W^'(fci?) 



a(k) 

Pl{1)in)ik,R) = bl.^W^kR). 



a2(fc) 



Pi2ik, R) is the sum of the following terms: 



Local 2P(1^o)(2o)(fc,i?) = 2 bio J Bo{k,q,-k- ti) 4'\-k~ q,q}W^{kR)W\qR) 



Local 2P('l^5)(2i)(fc,i?) - 61062 



10 



d\ BoiK q, -k - q_) _ ^ ^2(,^) 



(27r)3 
d3q 



,{k) 



(27r)3 



2P(1^)(2i)(^-^) = boih 
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2P{w)i22)ik,R) = hobn 
2^(n)(22)(fc,^) - boibn 
2PA^^.o,^(fc,i?) = 6infe02 



2Pii'lm){k,R) = boibo 



Bo{k, q, -k - q) W{kR) W{qR) W{\k + q\R) 
W{kR) W{qR) VK(|k + qji?) 
W{kR) W{qR) VK(|k + qji?) 
W{kR) W{qR) H/(|k + qJi?) 
W{kR) W{qR) VK(|k + qJi?) 
W{kR) W{qR) W^(|k + qji?). 



d3q 


5o(k,q, -k- 


q) 


(27r)3 


a(fc) 




d3q 


5o(k,q, -k- 


q) 


(27r)3 


a{q) 




d3q 


5o(k,q, -k- 


q) 


(27r)3 


a{k) a{q) 




d3q 


5o(k,q, -k- 


q) 



(27r)3 a(g)Q!(lk + q 
(fq Bo(k, q, -k-q 



(27r)3 a(/c)a(g)a(|k-t-ql) 



(77) 



(78) 



(79) 
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P22{k, R) is the sum of the following terms: 

Local, Gauss -P(2ci)(2o)(^' 



2 bio 



Local, Gauss 2P(2o)(2i)(^' -R) 



Local, Gauss P[-^){2i){k,R) 



2 6io &20 



(27r)3 
d3q 



(27r)3 
d3q 



Po(9) Po(|k - q|) 4^^(q, k - q) W{kR) W{qR) W{\k - q|ii) 



2 

2-P(20)(22)(^>-^) = 2610&II 



Po(<?) Po(|k - q|) v^2(|k _ q|i?) 

Po(|k - q|) 4^)(q, k - q) ^(^i?) W(|k - q|ii) 



(27r)3 

d'q Po(g) 



2-P(21)(22)(^'-^) 



(27r)3 

d3q Po(g) 



620611 / ^^Po{\^-c^\)W\qR)W\\k-<i\R) 



■phh 

-^(22)(22) 



Po(|k-q|) , Po((z) Po(|k-q|: 



.(|k-q|)2 a{q) a(|k-q|) 



H^2(gP) W^2(|k-q|P) 



2-^(20) (23) 



2-P(2i)(23)(^' — 620602 
2-F(22)(23)(^'-R) = 611602 



d3q Po(g) Po(|k- q|) , 9,, 



(27r)3 a{q) a(|k - q|) 

d3q ^^o(g)Po(|k-q|)^2(^^)^2(|k_q|^) 



(27r)3 a{q) a2(|k-q|) 



-^(23) (23) 



(80) 



Pi^{k, R) is the sum of the following terms: 



Local, Gauss 2P(^Q)(3o)(fc, P) 

2-P(n)(3o)(^'-^) 
Local, Gauss 2P(}o)(3i)(/;;, P) 

2-P('ll)(31)(^'-^) 
2-P(10)(32) (fc)-R) 
2-P(ll)(32)(^!-K) 
2-P(10)(33) (^)-K) 
2-P(ll)(33) (^)-R) 



2-P(io)(;m)(^>-^) 



= 6 6?oPc 



_d3q 
{2ny 



6 601 610 



Po(fc) 

a{k) 



Po{q)4''\Kci,-ci)W^{kR) 



(27r)3 

?3 



Po{q)4'\Kci,-ci)W\kR) 



46io62oPo(/c) y ^Po(9)ji'*^(q,k)W^(A;P)iy(5P)iy(|k + q|P) 



= 4601620 



Poik) 



a{k) 
2 6io6iiPo(A;) 

^'o(fc) 



4^ Po(g) J^^^ (q, k) W{kR) W{qR) W{\k + q|P) 
(Cf^ ^ '^^'^Cq, k) W{kR) WiqR) W(|k + q|P) 



2601 611 



a{k) 



4"' (q, k) W{kR) W{qR) W{\k + q|P) 



= 77610 621 Po (A;) 



(27r)3 aiq) 



(27r)3 a{q) 
d3q 



610 612 Po (A:) 
Poik) 



(2^)3 ^"('^^ 

2Pri^)(34)(fc,i?) = 3^01612^ / (03^OW 



1 



Q;2(g) 



+ 



a(/c) a(g) 
1 



W'^ikR) W^{qR) 
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Local, Gauss 2P^^^^^^,^^{k,R) = 610 630 ^oW j Poiq)W\kR)W\qR) 
'^Ptim5)ik,R) = boiho^ j -0^^Po{q)W\kR)W\qR) 

The only important two- loop contribution is (relevant for small /nl, large ^nl): 

P/'a^'^Cfc, R) = ^ 610630 j ^-^^ To(k, q, p, -k - q - p) WiqR) W{pR) W(|q + p|ii) - k - q - p|ii). (82) 

2. The halo-matter spectrum 

The full halo-matter cross spectrum at one loop is 

P''"'{k, z, R) = D^{z) P/r(A;, R) + D\z) P'^^(k, R) + D^{z) [Ps'Tl^, R) + Pi^ik, R)] • (83) 
Pii^{k, R) is the sum of the following terms: 

i (10)(10)^ 

Plwm){k,R) = boi^W^kR). (84) 



Local, Gauss Pp-^^^.^^Ak^R) = bw Po{k)W'^{kR) 



PvT'ik, R) is the sum of the following terms: 



Local 2P(1™)(20)(fc,^) = 2 6ioy ^ i3o(k, q, -k - q) j(^)(-k - q, q) ^^^(fcP) VF^^^^^ 
Ptm20)ik,R) = b,, J (0 ^"^^'aik^'""^ 4'\-).-^,ci)W\kR)W'{qR) 
Local P(1:5)(2i)(fc,P) = ^620 y^Po(k,q,-k-q)W^(A:P)l^(gP)l^(|k + q|P) 

ptw)i22)ik,R) = i&ii / ^(~)'' ~ ^(fe-^) ^(g-^) ^(i^^ + qi^) 

Ptm2S) ik, R) = lbo2j-^ X)<^(|k + q|f ^(^"^^ ^(^"^^ ^('"^ + ^1"^^- ^''^ 

-^22" (^5 i^) is the sum of the fohowing terms: 

/73 2 
^Po(g)Po(|k-q|) [j^^'(q,k-q)] W^{kR) 

Local, Gauss P^^^o) (21) {k,R) = 620 y ^ Po(g) Po(|k - q|) J^^)(q, k - q) W^(fcP) VF(gP) VF(|k - q|P) 

Pmi22)ik,R) = 6n / |^^Po(|k-q|)j(^)(q,k-q)W(fcP)W(9P)W(|k-q|P) 

Pt2m2S)ik,R) = bo2 J ^^^^^4'\ci,k-<DW{kR)W{qR)W{\k-ci\R). (86) 



Pi^{k, R) is the sum of the following terms: 



Local, Gauss 2P: 



<W^{kR) 



(11)(30) 



{k,R) = 3 601 



Po(fc) f d\ 



Po{q)4'\Kci,-ci)W'ikR) 



a{k) J (27r)3 



Local, Gauss P^{S)m{k,R) = 2b2oP{k) J Po{q) 4'\q, k)W{kR)W {qR)W{\k + q\R) 
Pmi32)(k,R) = 6nPo(fc) I -^^4'\ci,k)W{kR)W{qR)W{\k + ci\R) 
^(1^)(33)(fc,^) = lb2^Po{k) J ^^W^kR)W^qR) 



-P(10)(34)(^'-^) 



Local, Gauss -P(1™)(35)(A;, -R) 



= ^h2Poik) 



= IbsoPoik) J 



(27r)3 
(27r) 



Poiq) 



1 



1 



a^{q) a{k)a{q) 
^Po{q) W^kR) W^iqR) 



W'^{kR) W^{qR) 



(10)(36) 



{k,R) 



1, /• rf'q 

3 °^ a(fc) y (27r)3 a2(5) 



The only important two-loop contribution is (relevant for small /nl, large 5nl): 

PiT'"(fc, i?) = i 630 y To(k, q, p, -k - q - p) W{pR) VF(|q + p|i?) - k - q - p 



\R). 



